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We study the suitability of correlator product states for describing ground-state properties of 
two-dimensional spin models. Our ansatz for the many-body wave function takes the form of either 
plaquette or bond correlator product states and the energy is optimized by varying the correlators 
using Monte Carlo minimization. For the Ising model we find that plaquette correlators are best 
for estimating the energy while bond correlators capture the expected long-range correlations and 
critical behavior of the system more faithfully. For the antiferromagnetic Heisenberg model, how- 
ever, plaquettes outperform bond correlators at describing both local and long-range correlations 
because of the substantially larger number of local parameters they contain. These observations 
have quantitative implications for the application of correlator product states to other more com- 
plex systems, and they give important heuristic insights: in particular the necessity of carefully 
tailoring the choice of correlators to the system considered, and its interactions and symmetries. 



INTRODUCTION 



A. Background 



Modeling strongly correlated systems exactly for more 
than a few particles is not possible due to the rapid in- 
crease in the size of the Hilbert space with their number. 
However, there has been much success recently in using 
tensor network methods to numerically simulate strongly 
correlated systems P^ These approaches allow states with 
suitable properties, such as obeying an area law for their 
entanglement entropy, 2 to be efficiently represented by 
a network of tensors. They also provide efficient meth- 
ods of 'contracting' the network to allow expectation val- 
ues of operators, or more basically wave function ampli- 
tudes, to be calculated. In particular, matrix product 
states (MPS) can be used to accurately describe large sys- 
tems in one dimension and can be exactly and efficiently 
contracted.^ Crucially this means that expectation values 
can be calculated in a number of steps that only grows 
as a low degree polynomial with system size and ten- 
sor size. This has been crucial for the success of MPS 
based algorithms like the density-matrix renormalization 
group, 4 and time-evolving block decimationP Tensor net- 
works can be extended to higher dimensions using the 
projected entangled pair states (PEPS) construction^ 
but unlike MPS these do not permit efficient exact con- 
traction. Approximate contraction procedures in two di- 
mensions, such as those based on MPS methods f^ or 
tensor renormalisation group (TRG)p"^ must be em- 
ployed at the cost of having a much less favorable, though 
polynomial, scaling. Recently a procedure for perform- 
ing PEPS simulations combining TRG and Monte Carlo 
sampling has been proposed, and it shows promising re- 
sults for treating large bond dimensions.^ In contrast 
to the PEPS construction, both the multiscale entan- 
glement renormalisation ansatz (MERA)P^ and tensor 



tree networks (TTNs) 13 utilize a hierarchical structure 
where the tensors are interconnected by bonds according 
to a tree pattern, and can be contracted exactly and effi- 
ciently. Like PEPS, these methods scale as a high degree 
polynomial. 



An alternative less computationally expensive ap- 
proach to represent two-dimensional systems is to con- 
sider classes of states that while not being efficiently con- 
tract ible exactly, are efficiently and exactly samplable, 
i.e. for any given configuration the amplitude can be 
found exactly. The tensor network formalism provides a 
powerful framework from which to devise new and physi- 
cally tailored ansatze with these properties for which the 
closest approximate ground state can be found using vari- 
ational Monte Carlo methods. For example, string bond 
states use a product of overlapping MPS 'strings' and ac- 
quire their exact samplability from the contractibility of 
the underlying MPS. 14 An alternative ansatz is formed 
by building the wave function using superpositions of 
all possible coverings of singlets i.e. resonating valence 
bond (RVB) statefl^EJ, anc [ it has been shown that RVB 
states have a PEPS description. 19 Recently a new ansatz 
within this class has been proposed for simulating lattice 
systems: so-called correlat or p roduct states pIEIl or en _ 
tangled plaquette statespHHI which are equivalent and 
are hereafter referred to as CPS. A similar approach was 
also proposed in some earlier works ™ Correlator prod- 
uct states can be thought of as a more basic form of 
tensor network states, where correlations between sites 
are encoded explicitly in "correlator" building blocks so 
that the amplitude for each configuration is given by a 
product of their scalar elements. They provide a slightly 
simpler ansatz than string bond states, but with similar 
power and properties. 20 



B. Motivation 

The CPS construction has been applied to a variety 
of models, and in all cases the energies found compared 
well with those found using other methods, e.g. PEPS, 
MPS, stochastic series expansion (SSE), but with a much 
reduced computational cost. As we will discuss in detail 
below, the CPS description is incredibly versatile, and 
the correlator type can be varied without adversely af- 
fecting the complexity of the calculation. Thus, the type 
of correlator can be chosen to best match the properties 
or symmetries of the system. It is also possible to go be- 
yond PEPS, which are composed of short ranged bonds, 
and consider CPS that contain long ranged bonds and 
that are area-law violating. 2 

In previous work, the estimated ground state energy as 
a function of the plaquette size has been determined. 22 
What has so far been lacking is a systematic study of the 
relative behavior of different correlator types for describ- 
ing the important physical properties of different sys- 
tems. In particular bond correlators, where correlators 
are arranged to connect pairs of sites across the lattice, 
have not been studied in any detail. Bond correlators 
allow the possibility to directly encode long-range corre- 
lations between distant sites using a very small number 
of parameters. This is not possible with plaquette cor- 
relators, due to the exponential increase in the number 
of correlator elements with plaquette size, or with ma- 
trix product states, where distant correlations are medi- 
ated by intermediate nearest neighbor bonds. However, 
plaquette correlators do provide a more straightforward 
ansatz that imposes a less rigid short ranged structure. 
It is thus interesting to investigate what effect the choice 
of correlator product ansatz has on the ability to describe 
systems that possess long-range correlations, as well as 
the difficulties they may present during minimization. 

The CPS construction seems promising for describing a 
wide variety of spin models, and modeling these systems 
could provide insights for cold-atom simulations of quan- 
tum magnetism! 27 * 28 * In previous works it has also been 
noted that certain important states such as the Laugh- 
lin wave function} 2 ^ which cannot be efficiently described 
with other types of tensor network states, 30 can have an 
exact CPS description when bond correlators between all 
sites are included. 20 This illustrates that bond correla- 
tors can efficiently describe complex topological systems, 
and may be able to describe other fractional quantum 
Hall states p2 To aid future studies which will apply the 
CPS ansatz to more complicated systems, in this work 
we test their effectiveness using systems whose behavior 
is well-known. We examine in detail the performance of 
different CPS ansatze, but note that there is no guarantee 
that the state with the lowest estimate of the energy bet- 
ter reproduces any properties of the ground state, aside 
from of course energy, better than another with a higher 
energy estimate P^ Therefore, as well as estimating the 
ground state energy, we also investigate the critical be- 
havior, the long-range correlations and antiferromagnetic 



order, which provide useful benchmarks of their effective- 
ness. Our observations provide insight of quantitative 
and heuristic value into the accuracy and applicability of 
the CPS approach. 



C. Outline and main results 

In this paper, we investigate different correlator types, 
with the aim of identifying those that are most effective 
at capturing the behavior of the system, and find that 
the optimum correlator type is strongly dependent on 
the physical properties of the system. In section [TTJ we 
describe the different correlators types and outline the 
methods used to determine the ground state properties. 
In section Hill we examine the effectiveness of the different 
correlator types at describing the long range order and 
critical phenomena by applying them to the quantum 
transverse Ising model (TIM) on a square lattice. We 
find that bond correlators are better able to predict the 
critical point and represent expected long range correla- 
tions than plaquette correlators for the two-dimensional 
TIM, where bond correlators allow a larger range of sites 
to be covered for a given number of parameters and com- 



putational effort. In section IV we extend the treatment 
to different lattice geometries, finding similar behavior 
for the performance of the two correlator types. In sec- 
tion |V| we investigate the ability of plaquette and bond 
correlators to describe a more complex type of long range 
order and antiferromagnetism by applying them to the 
antiferromagnetic Heisenberg model (AFHM). We find 
that plaquette correlators are more successful at describ- 
ing the AFHM, since they provide a much larger number 
of local parameters to better capture the complicated lo- 
cal structure encoding the antiferromagnetic order in this 
system. Finally, we conclude and summarize the results 
in section IVTl 



II. METHOD 

The exact and efficient samplability of CPS makes 
them ideal for variational minimization (for details of 
the exact-efficient samplability of CPS see appendix [A| . 
Previously CPS have been minimized using a general- 
ized eigenvalue method,^ and a deterministic method. 21 
Here we use a Monte Carlo based stochastic minimiza- 
tion method to determine the ground stated which has 
previously been successfully applied to one-dimensional 
systems represented by matrix product states! 3 -^ The cor- 
relator elements that best approximate the ground state 
are found by estimating the derivative of the energy with 
respect to each correlator element, and then updating 
each correlator element according to the direction of the 
derivative with a random step size. This allows mini- 
mization using only the first derivative of the energy (so 
requiring fewer computation steps), and ameliorates the 
error associated with updating the correlator elements. 
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FIG. 1. (Color online) Illustration of the CPS ansatz for a 
state \ip) describing a simple 2x2 system with open boundary 
conditions. The local states Si are those of a spin-| system. 
Each correlator C l,J represents the amplitudes for different 
configurations of spins on two sites i,j. The four-site system 
has a total of sixteen basis states, with the weight for an 
example basis state W(UW = C\fC\fClfClf. 



Precise details of this numerical method are described in 
appendix [B] 

A. System set-up 

We consider a system with N sites where each site i 
has an identical local Hilbert space spanned by states \si). 
The full Hilbert space of the system is then spanned by 
states |s) = Isi)®!^)®* • " where we term s = (si, $2? * ' ' ) 
the configuration of the state. The many-body wave 
function is \ip) = ^ s Vl^(s)|s), where W(&) is the ampli- 
tude or weight of a given configuration s. In this work we 
will for simplicity consider only spin-^ systems, but the 
methods used can be applied straightforwardly to lattice 
systems possessing a larger on-site dimension. 

In the correlator product state description, the weight 
W(s) is given by the product of correlator elements Cs\ i} 
over the lattice 



W(s) 






(1) 



where each correlator element is a c-number that de- 
scribes the amplitude of a configuration s^y of a sub- 
group of sites {i}. The wave function in the CPS repre- 
sentation is given by 



iv>) = EII^wi s i---^>- 



(2) 
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FIG. 2. (Color online) Illustration of the different correlator 
types used. The blue circles represent lattice sites, and only 
correlators that link these 3x3 sites are represented, (a) 
Nearest neighbor (n. n.) bond correlators, (b) Bond correla- 



tors with r m 

2x2 plaquette correlators. 



1. (c) Bond correlators with r n 



2. (d) 



We illustrate this in Fig. [T] for the case where the sub- 
group of sites is a nearest neighbor pair i.e. 'bond' corre- 
lators, for a 2 x 2 system. 

This compact description of a state can be extended to 
any form of correlator, for example each correlator could 
represent four sites in a plaquette, or a string of a given 
number of sites. One advantage of using CPS is that 
the wave function amplitude is a simple product of c- 
numbers. Unlike with PEPS, the amplitude W(s) can be 



calculated efficiently and exactly for any configuration s. 
This allows the wave function to be efficiently sampled, 
so that Monte Carlo methods can be used to determine 
expectation values and minimize the wave function. The 
flexibility and exact-efficient samplability of CPS is de- 
scribed further in appendix [Aj 



B. Correlator types 

In the following, we study the ground states using dif- 
ferent types of correlator, and compare their properties. 
We use periodic boundary conditions with the transla- 
tional invariance of the system allowing one correlator of 
each type to describe the entire system. 

The types of correlators used for square lattices are 
illustrated in Fig. |2J Nearest-neighbor (n. n.) bond cor- 
relators, with one correlator for the vertical bonds and 
one correlator for the horizontal bonds as shown in Fig. 
|2|a), are the simplest correlator types used. This gives a 
description of the ground state using only eight param- 
eters. We also use bond correlators with longer range 
bonds. In general we term a bond correlator of size 
r max as including all bond vectors up to the maximum 
range Ar max = (r max ,r max ) i.e. it includes all the diag- 
onal bonds as well as bonds along the lattice axes. Fig- 
ures f2^b) and (c) show the set-up with correlators of size 
r max = 1 and r max = 2 respectively, and in the following 



calculations we use correlators up to r max = 8. For cor- 
relators of size r max , the number of correlator elements is 
8r max (r max + 1). To calculate the energy, the number of 
computational steps scales as 0(A/" 2 r^ ax ), and the only 
part of the calculation that depends on the correlator size 
and type is the calculation the correlator fraction given 
in equation (B2). 



We also use plaquette correlators, with the smallest 
2x2 plaquette illustrated in Fig. [2^d), up to a size of 
4x4 plaquettes. The correlators are set up so that they 
are displaced from one another by one site and overlap 
one another. In general, the greater the overlap between 
the correlators the more accurate the description of the 
ground state. The number of elements in an n x n cor- 
relator scales as 2 n , so the memory requirements for a 
given plaquette correlator size scale much faster than for 
bond correlators. However the number of computational 
steps required to calculate the energy is 0(7V 2 n 4 ), which 
is a mild polynomial scaling and grows much more slowly 
than the number of correlator elements. The reason for 
this is that calculating an off-diagonal expectation value 
only involves picking out the correct correlator element 
for the subgroup of sites spanned by a correlator. As a 
result it depends on the number of sites spanned by the 
correlator and the number of correlators that fall on a 
given site, and not on the number of elements. Note that 
calculation of the energy derivative formally requires a 
number of steps that does scale with the number of cor- 
relator elements, however in practice this part of the cal- 
culation is fast compared to estimating the energy itself 
and it is not rate limiting. Thus plaquette correlators al- 
low the use of many more parameters for describing the 
system with only a modest increase in computational ef- 
fort, although in our calculations their size is ultimately 
limited to 4 x 4 by memory requirements. 

Comparing the computation time for correlator types, 
a calculation using 3x3 plaquette correlators has ap- 
proximately the same number of computation steps as for 
bond correlators up to r max = 4, and a calculation using 
4x4 plaquette correlators has approximately the same 
number of computation steps as for bond correlators up 
to r max = 8. However, when comparing the performance 
of different correlator types, it is worth considering the 
three main differences between them: (i) The number 
of sites that can be reached from a given site using the 
correlator: r max = 1 is equivalent to a 2 x 2 plaquette, 
r max = 2 is equivalent to a 3 x 3 plaquette, and r max = 3 
is equivalent to a 4 x 4 plaquette. (ii) The computational 
effort: for a given number of spanned sites, our routines 
for bond correlators are 2, 3.5 and 5.5 times faster than 
for plaquette correlators respectively, (iii) The number 
of correlator elements: the number of correlator elements 
for a given number of sites spanned is far smaller for bond 
correlators. In addition to the reduction in computa- 
tional effort, this makes it possible to span far more sites 
using bond correlators. The more fragmented structure 
of bond correlators also permits different minimization 
strategies as described in more detail in appendix [Bl 



III. LONG RANGE CORRELATIONS AND 

CRITICAL BEHAVIOR: THE QUANTUM 

TRANSVERSE ISING MODEL ON A SQUARE 

LATTICE 

As an ideal test case for examining different correlator 
types, we consider the TIM on an L x L square lattice, 
described by the Hamiltonian 



H = -J 



i=i,j=i 



( 4«v 



»+i,j] 



■^Vp 1 ^^), (3) 



where J > is the coupling, i and j denote the lat- 
tice index in the two perpendicular directions, cr?^ , is 
the Pauli z(x,y) operator for spin (i, j), and g is the di- 
mensionless transverse magnetic field. We consider this 
model because it is one of the archetypal systems that 
exhibits a quantum phase transition at a finite magnetic 
field. 34 The system moves from an ordered ferromagnet 
in the z direction at g = to a state disordered in the 
z direction and aligned in the x direction when g ^> 1. 
As the system approaches criticality at g = g c , the corre- 
lation length diverges, and the system becomes gapless, 
making a numerical description of the state challenging. 
The critical point in the thermodynamic limit, found us- 
ing a finite-size scaling analysis, 35 is g c « 3.044, where 
the pseudo-critical point 36 was defined using a careful 
extrapolation of the ratio of the energy gap between con- 
secutive system sizes. 37 

By applying the CPS ansatz to this system, we inves- 
tigate its performance at describing the behavior of sev- 
eral important physical quantities as the transverse mag- 
netic field is varied. We focus in particular at modeling 
the long-range correlations close to the pseudo-critical 
point, and compare with the many previous studies per- 
formed using other numerical and (approximate) analyt- 
ical methods. This system is thus a highly useful bench- 
mark of the method for demonstrating both the effective- 
ness and limitations of the CPS approach. 



A. Energy 

The energy and its derivative, as well as other observ- 
ables, are calculated using the algorithm described in ap- 
pendix [B| Note that for all correlator types we restrict to 
only real parameters without any loss of generality, since 
for this model the ground state can be constructed using 
real, positive weights for all configurations. We first cal- 
culate the energy at a value of the magnetic field close 
to the critical point for two systems. The largest sys- 
tem that has been solved numerically exactly is a 6 x 6 
system,^ and so provides a good comparison for how 
well the method is working. This system size has also 
been solved using TTNsp^ so we also compare the ac- 
curacy of the CPS method to those results. The en- 
ergy of the ground state found using different correla- 
tor types is shown in Table [TJ The error, defined as the 



Correlator type 


No. of elements 


Energy 


Error 


n. n. bonds 


8 


-3.2348(4) 


13 x 10" 3 


'"max — -L 


16 


-3.2384(3) 


9 x 10" 3 


^*max == -^ 


48 


-3.2457(3) 


2 x 10" 3 


''max — O 


96 


-3.2465(2) 


1 x 10" 3 


2x2 plaquettes 


16 


-3.2387(4) 


9 x 1(T 3 


3x3 plaquettes 


512 


-3.2463(2) 


1 x 10" 3 


4x4 plaquettes 


65,536 


-3.24725(5) 


2 x 10" 5 



TABLE I. Energy per site in units of J found using stochastic 
minimization for the two-dimensional TIM in a 6 x 6 system 
for different correlator types. Results are compared with ex- 
act results at g = 3.05266, which has been calculated to be 
the pseudo-critical point in a 6 x 6 system, with a ground 
state energy of -3.2472744 • • • (Ref.l55t. 



Correlator type 


No. of elements 


Energy 


n. n. bonds 


8 


-3.2325(1) 


'"max — -L 


16 


-3.2357(1) 


''max — -^ 


48 


-3.23815(8) 


'"max — O 


96 


-3.23865(5) 


'"max — 


240 


-3.23866(7) 


'"max — O 


576 


-3.23882(7) 


2x2 plaquettes 


16 


-3.23599(9) 


3x3 plaquettes 


512 


-3.23864(4) 


4x4 plaquettes 


65,536 


-3.23903(5) 



TABLE II. Energy per site in units of J found using stochastic 
minimization for the two-dimensional TIM in a 31 x 31 system 
for different correlator types at g — 3.05. 



difference between the energy estimated using CPS and 
the exact ground state energy, is also displayed. As ex- 
pected, we find that the energy converges to the exact 
value for increasing correlator size. In comparison, for 
TTNs, around 10 7 parameters are required to reduce the 
error to 2 x 10 -5 . 13 We also apply the CPS method to 
much larger systems. In Table [TT| we show the minimized 
energy in a 31 x 31 system at g = 3.05, which is close to 
the pseudo-critical point. 

We find that even though plaquette correlators have 
more elements, the energy convergence is more well- 
behaved than for bond correlators. In line with com- 
monly known properties of non-convex optimization, it is 
likely that having a larger number of parameters is ben- 
eficial when there is a foliated energy landscape since it 
allows the optimization more freedom in finding the min- 
imum energy state. Conversely constraining the number 
of parameters for such long range bonds tends to result in 
numerous local minima giving a more difficult minimiza- 
tion problem. For example in a previous work minimiz- 
ing over a similar class of many-body states, it was found 
that the enforcement of certain symmetries increased the 
difficulty of finding a good state, since this amounted 




FIG. 3. (Color online) Difference in energy between ground 
state calculated using nearest-neighbor bond correlators, and 
ground state calculated using 4x4 plaquette correlators. The 
dashed line indicates the critical magnetic field found by a 
finite-size scaling analysis (Ref. 135)) . and the dotted line in- 
dicates the peak of the energy difference at g = 3.17 — this 
is approximately midway between the pseudo-critical points 
found with nearest-neighbor bond and 4x4 plaquette corre- 
lators. 



to cutting through the energy landscape of the parame- 
ter space, dividing it into separated minima. 38 However, 
even though a larger span of sites is found to be needed 
to minimize to a given energy with bond correlators, the 
number of parameters required is still far smaller than 
for plaquettes substantially reducing the computational 
effort. 

We also investigate how the choice of correlator affects 
the energy at different values of the transverse magnetic 
field g. Figure [3] shows the difference between the energy 
calculated using n. n. bond correlators (which give the 
highest estimated ground state energy) and the energy 
calculated using 4x4 plaquette correlators (which give 
the lowest estimated ground state energy), for different 
values of the transverse magnetic field g. There is not 
much difference between the two values for \g — g c \ ^> 1, 
but the energy difference increases when g « g c . Far 
from criticality, when the correlations are expected to 
be short range, the energy can be calculated accurately 
using a small number of parameters, however as the cor- 
relation length increases a larger number of parameters 
are needed to describe the system accurately. The max- 
imum energy difference is still small (< 1%), however 
this can lead to large differences in the properties of the 
state due to the large number of low-lying excited states 
combined with a vanishing gap as g c is approached. This 
is a general problem with using a variational approach 
to describe critical systems. The investigation of these 
properties is described in the following sections, and we 
find that choosing an ansatz with a suitable structure 
helps to better describe some important physical prop- 
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FIG. 4. (Color online) Absolute magnetization (\cr z \) as a 
function of transverse magnetic field g in an L x L system, 
for the ground state found using 2x2 plaquette correlators. 
The pseudo-critical point occurs at g w 3.15. Inset: (\a z \) 
at constant g > g c as a function of L, both axes having a 
logarithmic scale, with a fit of the form (\<j z \) oc L~ b . The 
data was fitted with b — 1.0. 
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FIG. 5. (Color online) Absolute magnetization [(a) and (b)] 
and transverse magnetization [(c) and (d)] as a function of 
transverse magnetic field in a 31 x 31 system, for the ground 
state found using different sized plaquette correlators [(a) and 
(c)] and bond correlators [(b) and (d)]. The dashed line indi- 
cates the critical magnetic field found by a finite-size scaling 
analysis in Ref. 35, which occurs at g c = 3.044. 



erties of the ground state even close to criticality. Even 
though formally the energy found using bond correlators 
is larger the structure of the ansatz seems to favor those 
states that possess long-range correlations. 



B. Transition point 

We investigate the position of the pseudo-critical point 
in the system by examining the local order parameters. 
Specifically, we calculate both the transverse magnetiza- 
tion (a x ), and the absolute magnetization (|cr 2 |), defined 
as the expectation of the absolute value of the average of 
all L x L spins, i.e. it quantifies how well the spins are 
aligned with one another (N.B. (a z ) is zero due to the 
global Z2 symmetry of the system) . Figure [4] shows the 
results for the absolute magnetization (\cr z \) as a function 
of transverse magnetic field using 2x2 plaquette corre- 
lators for a number of L x L system sizes. The error bars 
for each point, where the error is given as the standard 
deviation of the G different bins (see appendix [B| , are 
smaller than the marker for the data point. 

As expected, we find a sharp drop in the magnetiza- 
tion at around g = 3.05 and the change in magnetization 
becomes steeper as the system size increases. We also 
find that the small non-zero magnetization for g > g c de- 
creases as the system size increases: for L = 51, the mag- 
netization at g = 4 is (\cr z \) w 0.022. The inset in Fig. [4] 
shows the magnetization as a function of lattice size for 
different fixed values of g > g ci with a fit to a power law 
decay of the magnetization, i.e. (\cr z \} g>gc = aL~ b . The 



data fit an inverse scaling with system size, indicating 
that the magnetization decays to zero for g > g c in an 
infinite system. While the behavior found for 2 x 2 cor- 
relators is not quantitatively precise (the pseudo-critical 
point occurs at g ~ 3.15 rather than at g = 3.044), it is 
simple and offers a computationally affordable means of 
approximately locating a critical point once the appro- 
priate local order parameters signifying it are known. 

As described in section [Til A[ we find that larger cor- 
relators are better able to minimize the energy. Figure 
5 shows how the behavior of the local order parameter 
(\<J Z \) and the transverse magnetization (a x ) depend on 
the correlator type and size. The position of g c , the crit- 
ical point found by a finite-size scaling analysis,^ is also 
indicated. The plots of the absolute magnetization (\cr z \) 
shown in Figures |5|a) and (b) show that as the correlator 
size increases, the calculated pseudo-critical point moves 
closer to g c for both plaquette and bond correlators. The 
transverse magnetization (shown in Figures [5Yc) and (d)) 
displays the expected 'knee' at the pseudo-critical point, 
which again moves closer to g c as the correlator size in- 
creases. The magnetization close to g c has not completely 
converged for the larger plaquette correlators in our cal- 
culations: the difference in (a x ) for the 4x4 plaquette 
correlator and the 3x3 plaquette correlator is 0.003. 
Comparing this with the convergence found using TTNs, 
for a 10 x 10 system, the difference in the transverse mag- 
netization for ~ 10 6 parameters and ~ 10 T parameters 
is around 0.002™ We see some scatter in the behavior 
of both the absolute and transverse magnetization calcu- 
lated using bond correlators, owing to the more difficult 
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FIG. 6. (Color online) The connected two-point correlation 
function plotted using a logarithmic color scale in a 31 x 31 
system at g = 3.05 found using (a) n. n. correlators and (b) 
bond correlators with r ma x = 8. (c) Correlation length as a 
function of g for different correlator types, bond correlators 
are open circles, plaquette correlators are closed circles. See 
Fig. [5] for full legends. The dashed line indicates the critical 
magnetic field found by a finite-size scaling analysis (Ref . l35|) . 



minimization, which means that the convergence cannot 
be accurately determined. However, despite this noise it 
is clear from our results that bond correlators display a 
pseudo-critical point much closer to that found in Ref. [35] 
for r max > 5 as indicated by both the absolute and trans- 
verse magnetization. This suggests that bond correlators 
are capturing the critical behavior of the local order pa- 
rameters in the TIM better. It is perhaps surprising that 
the energy density is better described by plaquette cor- 
relators, while the behavior of the local order parameter 
is qualitatively better with bond correlators. We next 
study the performance of plaquette and bond correlators 
at describing the long-range correlations, to provide some 
further insights. 



C. Long-range correlations 

The success of the MPS and PEPS-type tensor network 
approach is intimately connected with the decay of cor- 
relations in the system. When the correlations are short 



range, the tensor network representation is likely to be 
able to model the system well, with a number of param- 
eters exponentially smaller than the Hilbert space size. 
When the long-range correlations approach polynomial 
decay, however, it can be difficult to model the system 
accurately using these methods, and a larger number of 
parameters are required to describe the entanglement in 
the system. Conversely, the success of MERA and TTNs 
is based on the logarithmic rather than linear scaling of 
the distance i n the netw ork between two points a given 
distance apart [ffl 13 l 39 l 4Q l anc [ this ability to describe long 
range correlations can dramatically aid performance. 

To investigate how the success of the CPS description 
is related to the ability to describe long range correlations 
we calculate the connected two point correlation function 
given by 



C x , y — (cr. 
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■ [x,y] 



I), (4) 



for a wave function that has been minimized at different 
values of the transverse magnetic field and for different 
correlator types. The connected two point correlation 
function C XjV describes the probability of two spins on 
separated sites being aligned with one another, and can 
only be non-zero if entanglement exists between these two 
sites in the underlying ground state. When the system 
is close to the pseudo-critical point, it possesses its max- 
imum correlation length £, and the system is the hardest 
to simulate numerically. 

We calculate the average connected two-point corre- 
lation function in a 31 x 31 system. Figures [6la) and 
p[b) show a plot of C x , y for n. n. bonds and bonds 
with r max = 8 respectively. These show that the cor- 
relations grow substantially with only a moderate bond 
length increase. We also determine the correlation func- 
tion Ci = C XjV as a function of distance I between the 
central point (0,0) and the point (x,y) and calculate 
the correlation length £ by fitting to an exponential de- 
cay C\ = exp (—//£) for different correlator types. We 
plot this as a function of transverse magnetic field g in 
Fig. |6jc) and the behavior illustrates the longer reach of 
bonds. We see a peak in the correlation length, consis- 
tent with the transition point seen in the magnetization. 
This peak becomes sharper and moves from g ~ 3.15 
for smaller bond correlators and plaquette correlators 
to g ~ 3.05 for the largest bond correlators, which is 
closer to the predicted value in the thermodynamic limit. 
However, even with the largest bond size, the correlation 
length does not exceed significantly more than a single 
lattice site, and thus seem to poorly model the true di- 
vergence in this property. It is clear that using nearest- 
neighbor bonds to describe the system discards much of 
the information about long-range correlations in the sys- 
tem. The increasing size of the peak gives some indica- 
tion of a slowly growing divergence, with bond correlators 
having r max = 8 showing the largest correlation length. 

The above results show that bond correlators are not 
only better for capturing the signatures of critical behav- 
ior as reflected in the local order parameters, they are also 




FIG. 7. (Color online) Correlators used in triangular and 
hexagonal lattice geometries, (a) Four-site plaquette correla- 
tors on a triangular lattice, (b) n. n. bond correlators on a 
triangular lattice, (c) r ma x = 1 bond correlators on a triangu- 
lar lattice, (d) Six-site plaquette correlators on a hexagonal 
lattice. 



better at capturing the long-range correlations despite 
having a larger energy estimate than plaquettes. Note 
that it is not uncommon for a state with a higher energy 
to represent other properties of the ground state more 
accurately than that with the lower energy estimated 
For this reason it is important that the choice of ansatz 
should try to reflect some expected underlying properties 
of the ground state and that the results, beyond just en- 
ergy, need to be carefully examined (as in our work). The 
energy of a given state depends entirely on short-range 
correlations. Plaquette correlators have a large number 
of parameters to describe the short-range interactions, 
while the comparatively small number of parameters for 
bond correlators limit the accuracy to which it can de- 
scribe the energy. However the structure of bond cor- 
relators allows an efficient description of any prevalent 
long-range correlations since it provides direct bonds be- 
tween more distant lattice sites than can be reached with 
plaquette correlators. Although the calculated correla- 
tion length is still not more than two lattice sites close to 
the pseudo-critical point, it is for this reason that bond 
correlators are better able to describe the longer ranged 
properties of the Ising ground state, as found in this in- 
vestigation. 



IV. DIFFERENT GEOMETRIES: ISING MODEL 

ON TRIANGULAR AND HEXAGONAL 

LATTICES 

To generalize the above conclusions to different lattice 
geometries, we also apply the CPS ansatz and minimiza- 
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FIG. 8. (Color online) Absolute magnetization as a function 
of transverse magnetic field for lattices with different num- 
bers of sites TV in (a) a triangular lattice with the ground 
state found using four-site plaquette correlators and (b) a 
hexagonal lattice with the ground state found using six-site 
plaquette correlators. The dotted line in both figures shows 
the position of the critical point calculated using a finite-size 
scaling analysis (Ref. 141 j) . 



tion method to triangular and hexagonal lattices using 
the correlator types illustrated in Fig. [7J The behavior 
of the system is expected to be qualitatively the same for 
the square lattice, with the critical magnetic field shifted 
due to the different coordination numbers. For hexag- 
onal and triangular lattices, finite-size scaling analysis 
predicts g c = 2.13 and g c = 4.77 respectively.^ We cal- 
culate the absolute magnetization as a function of trans- 
verse magnetic field for the smaller plaquette correlators, 
for different lattice sizes, and the results are shown in 
Fig. [8j We see the expected behavior, with the pseudo- 
critical point shifted to a higher magnetic field for small 
plaquette correlators, and find that the value of the ab- 
solute magnetization at g ^> g c decreases for increasing 
lattice size L. 

As with the square lattice, we calculate energy close 
to the pseudo-critical point, and we find that plaquette 
correlators give the lowest energy. For example, for a 
triangular lattice with 400 sites at g = 4.77 the energy 
found using 4x4 plaquettes is —4.9899(1) J, while with 
bonds having r max = 5 the minimum energy is found 
to be —4.9895(2) J. Also as in square lattices we find 
that despite this, bond correlators seem to be better for 
determining the critical behavior of the system. Figure 
|9Fa) shows the absolute magnetization as a function of 
transverse magnetic field for different correlator types. 
As with the square lattice, the longer range of the bond 
correlators predict a pseudo-critical point that is closer 
to the critical point (although once again more scatter is 
evident). The correlation length as a function of trans- 
verse magnetic field in shown in Fig.^b). Again, we find 
that the largest correlation length for bond correlators, 
and the position of the pseudo-critical point is consistent 
with that predicted by the local order parameter. 

The numerical calculations indicate that the results 
found for the TIM are not dependent on the lattice ge- 
ometry chosen. For all the geometries investigated it is 
found that the longer ranged bond correlators predict a 
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FIG. 9. (Color online) (a) Absolute magnetization as a func- 
tion of transverse magnetic field in a 400-site triangular lat- 
tice, for the ground state found using different- sized plaquette 
correlators (closed circles) and bond correlators (open circles), 
(b) Correlation length as a function of transverse magnetic 
field found using different-sized plaquette correlators (closed 
circles) and bond correlators (open circles) . The dotted line in 
both figures shows the position of the critical point calculated 
using a finite-size scaling analysis fRef. I4lj) . 



critical point closer to that calculated in previous works, 
and can better describe the long-range correlations that 
are expected close to criticality, despite the lower mini- 
mized energy found using plaquette correlators. 



V. ANTIFERROMAGNETIC ORDER: THE 
HEISENBERG MODEL 



To examine how plaquette and bond correlators com- 
pare when describing very different physics, we also in- 
vestigate the AFHM on a square lattice. This system is 
more challenging to describe than the TIM, as it exhibits 
richer behavior owing to its symmetry and antiferromag- 
netism. The antiferromagnetic spin- 1 , Heisenberg model 



is described by the Hamiltonian 

L 

H = J J2 S^-S^ +1 '^H-S^].S^ +1 ], (5) 

where S = \((J X1 (J yi cr z ) is the spin-^ operator. It cor- 
responds to the half-filled band limit of the Hubbard 
model, and is assumed to describe the antiferromagnetic 
undoped insulator La2Cu04 as well as other undoped 
copper oxide materials.^ 

For a pair of spins, the ground state of Ho is simply 
a spin singlet. However, due to the monogamy of en- 
tanglement, each spin cannot be in a singlet state with 
more than one neighbor, and the ground state for N spins 
has to satisfy both local minimization and global trans- 
lational symmetries giving a highly complicated super- 
position stateP As a result of this, it is predicted that 
long-range order arises!™ This is characterized by the 
two-point correlation function for the total spin (S^ -S^) 
and the staggered magnetization, which is reduced from 
the Neel, or classical, value by quantum fluctuations P^ 

The AFHM has also already been treated using plaque- 
tte type correlators p° l 22 l but here we extend that descrip- 
tion by comparing them with bond correlators, and ex- 
amine which correlator type better captures the expected 
long-range order. The results obtained here are also com- 
pared with those found using SSE, which give the cur- 
rent best estimate. 44 This system has been treated using 
a resonating- valence bond (RVB) picture, by considering 
states that are superpositions of all possible coverings of 
singlets on the sitesP^"^ The latest study 18 gives esti- 
mates of the energy and staggered magnetization that 
agree well with those found using SSE. 

We restrict to configurations with an equal number 
of up and down spins, and move between configura- 
tions by flipping a pair of opposing spins (so that the 
net magnetization is unchanged). This can be done 
since the ground state exists in the spin sector having 
J2 { ■ (Jz = 0.™ We assume real, positive correlators 
and use the Marshall-sign rule to determine the correct 
sign for each configuration. 45 Otherwise we apply the 
same minimization routine as used above for the TIM. 
The ground state energy was found for different plaque- 
tte types and lattice sizes, and the results for a 14 x 14 
lattice are shown in Table IIIII 

For plaquette correlators, the energies obtained are 
consistent with previous calculations! 20 1 22 1 We see that 
plaquette correlators are able to minimize the energy 
much better than bond correlators, and that the discrep- 
ancy between the correlator types is larger than with the 
TIM. For example, for L = 14, even 3x3 correlators 
are better able to minimize the ground state energy than 
bonds with r max = 7, where bonds are long enough to 
couple every pair of sites on the lattice, and where the 
number of parameters is comparable. Thus the relative 
behavior of bond correlators with respect to plaquette 
correlators is worse than with the TIM. This is likely to 
be as a result of a more complicated energy landscape for 
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Correlator type 


No. of elements. 


Energy 


^max — -L 
^max — ^ 
'"max — O 
^max — ^ 
^max — O 
'"max = t> 
''max — < 


16 
48 
96 
160 
240 
336 
504 


-0.6560(3) 
-0.6627(3) 
-0.6641(2) 
-0.6646(1) 
-0.6648(3) 
-0.6647(2) 
-0.6649(3) 


2x2 plaquettes 
3x3 plaquettes 
4x4 plaquettes 


16 

512 

65,536 


-0.6567(3) 
-0.6662(2) 
-0.6685(4) 



TABLE III. Energy per site E/J found using stochastic min- 
imization for the AFHM in a 14 x 14 system for different 
correlator types. The result from using a SSE is E/J — 
-0.670222(7) (Ref.[44j). 
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FIG. 10. (Color online) Staggered magnetization found us- 
ing different plaquette types against the inverse of the lat- 
tice size L. (a) shows the staggered magnetization M? = 
Cl/2,l/2 and (b) shows the staggered magnetization Mf (L) = 
jj2 J2i j X~ iy +J Ci,j- I n both plots the black dashed line 
denotes the staggered magnetization in the thermodynamic 
limit found using SSE fRef.l44l). 



the AFHM than for the TIM, such that the fragmented 
structure imposed by the bond correlators makes it dif- 
ficult to describe the ground state. Consistent with the 
RVB picture, it therefore appears that having a large 
number of local parameters is more important to accu- 
rately describe the ground state of this model than to 
provide longer-ranged parameters. This is further high- 
lighted by the fact that the energy estimated using bond 
correlators with r max = 4 is the same (within error) as 
the energy estimated with r max = 7, so the improvement 
with increasing bond size is seen to saturate. 

As with the TIM, we also investigate how well the cor- 
relations are described by the different correlator types 
by calculating the staggered magnetization. The stag- 
gered magnetization has been previous ly ca lculated us- 
ing a variety of tensor network methods I^JiMMI including 
correlator product statesP^It has been found to be chal- 
lenging to determine accurately, with the estimated value 
in general being larger than the value found using SSE. 44 
The staggered magnetization M\(L) can be defined using 
the two-point correlation function at the furthest point 
out in the lattice: 49 



M 2 1 (L)=C L/2 , L/2 = {S^-S^ 2 ' L ^}. 



(6) 



An alternative definition for the staggered magnetization 
ED 



is: 



M\{L) 



i£(-ir + ^ 



L 2 



(7) 



These scale to the (same) staggered magnetization in the 
thermodynamic limit M(oo). The scaling behavior of 
Mi ^ has been predicted using chiral perturbation the- 
ory and spin- wave theory, and is given by M\ 2 (L) = 
M(oo) 2 + 6 1 , 2 /L + --..E 4 HniIIl 

We calculate the two-point correlation function, and 
similar to a previous workf^ it is found that the z- 
component of the correlations departs from the x, y- 
component of the correlations — the z correlations tend 



to zero while the x, y correlations remain finite. This 
is likely to be due to the preferential treatment of the 
z axis in the Monte Carlo method and inherent in the 
CPS correlator elements, i.e. the chosen basis. We also 
find that very small changes in the estimated ground 
state energy can lead to large differences in the estimated 
staggered magnetization. The results of these calcula- 
tions are shown in Fig. 



10 In contrast to the results 



in previous sections where bond correlators performed 
better, we find that plaquette correlators provide a bet- 
ter estimate of the long-range order, with an extrapo- 
lated value of Mi(oo) = 0.33(1) and M 2 (oo) = 0.322(7) 
found using 4x4 plaquette correlators closest to the 
value M(oo) = 0.3070(3) found using stochastic series 
expansion. 44 The value found in a previous study using 
the RVB basis is M(oo) = 0. 30743(1) P Unlike with the 
TIM, we find that it is not the range of the correlators 
that determines the accuracy with which correlations can 
be represented, but instead seems to be the number of lo- 
cal parameters. Note that although the RVB model has 
a similar pairwise structure to bond correlators, the de- 
scription used here is not equivalent to the RVB model, 
which does seem better able to represent the ground 
state. This is most likely due to the singlet structure 
of the RVB model, which can preserve the symmetry be- 
tween the three spin axes exactly. The relative behavior 
of bond and plaquette correlators is consistent with the 
minimized energy, and highlights the importance in the 
Heisenberg model of describing accurately the more com- 
plicated local correlations between multiple nearby spins 
in the three spin axes. 



VI. SUMMARY AND OUTLOOK 

We have investigated the performance of the most com- 
mon forms of correlator product states for describing the 
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ground state of the TIM in various two-dimensional ge- 
ometries. Even correlators with a small number of pa- 
rameters (16 parameters for the 2x2 plaquette correla- 
tor) can describe the qualitative behavior of large systems 
(up to 51 x 51 spins). Increasing the correlator size to 
4x4 site plaquettes or bonds with r max > 3 provides 
accurate results for up to L = 31 (and in principle larger 
systems could be studied). We have also found that, 
despite having a slightly larger minimized energy, bond 
correlators are more successful than plaquette correlators 
at describing critical behavior, such as the long-ranged 
correlations, whilst using a smaller number of param- 
eters and having a faster computation time. We have 
also investigated the AFHM in a square lattice. Results 
with plaquette correlators are comparable with previous 
works ™ and we have found that, unlike the TIM, pla- 
quette correlators are better able to describe both the 
energy and the long range behavior of the ground state. 

The system sizes reached here are comparable with 
experimental cold-atom lattice systems, and the abil- 
ity to efficiently model spin systems is useful for 
rapidly progressing cold-atom simulations of quantum 
magnetism! 27 * 28 * A general result for the optimal corre- 
lator structure for a given system is not straightforward 
to deduce, but our findings still provide useful insights for 
applications to condensed matter systems. For example, 
we would expect plaquette correlators to be the optimal 
choice for the J\ — J2 model studied in Ref. [22] and 
for the Fermi-Hubbard model. However, we may expect 
bond correlators to perform better for the Bose-Hubbard 
model, where the choice of the local Fock basis may al- 
low a good description of the ground state with a smaller 
number of local parameters, or in spin-systems which are 
highly anisotropic like the TIM. In addition, the Laugh- 
lin wave function^ has an exact CPS description using 
bond correlators. 20 

It may be that an optimal CPS ansatz consists of a 
hybrid of bond and plaquette correlators, such as strings 
of sites (in a similar set up to string bond states) linking 
all sites on the lattice, or extended rectangular plaquette 
correlators. This would provide a larger number of local 
parameters while still allowing the correlators to have a 
long range. In further work, 52 the CPS description will 
be applied to other systems, such as the Bose-Hubbard 
model in the presence of an artificial magnetic fieldp^ in 
the regime where the fractional quantum Hall states are 
predicted to occur! 5 -^ 



Appendix A: Exact and efficiently samplable tensor 
network states 



For the case of n. n. bond correlators, the resulting CPS 
is in fact equivalent to an MPS with internal dimension 
X equal to the physical dimension dW* Note that this 
equivalence does not hold in the reverse direction, i.e. 
not all MPS with x — d are equivalent to a CPS. Corre- 
lator product states can be thought of as a special class 
of matrix product states that can be factorized using the 
copy tensor (also known as the diagonal in category the- 
ory, the COPY-gate or the COPY-dot in circuits)!^ The 
copy tensor copies a chosen set of physical basis states 
e.g. \x) where x =t,i, and subsequently breaks up into 
disconnected states. This is essential for such states to 
be efficiently and exactly sampled. 
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FIG. 11. (Color online) Illustration of the equivalence be- 
tween matrix product states and correlator product states for 
the simple case of bond correlators, when the dimension of the 
correlators is the same as the physical dimension: (a) Shows 
the correlators connected by copy tensors (black dots), which 
is equivalent to an MPS where (i) each MPS tensor (dashed 
ellipse) is formed by the copy tensor and the neighboring cor- 
relator or (ii) each correlator is 'split' (e.g. using a singular 
value decomposition) and each MPS tensor (dashed box) is 
formed from a copy tensor and the adjacent correlator fac- 
tors; (b) this can be factorized, breaking up the connections 
while ensuring the same configuration is present on neighbor- 
ing correlators; (c) this gives a product of scalar correlator 
elements. 
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This is illustrated in Fig.[TTJa), where the copy tensor 
is represented using a black dot, while the large circle 
represents a matrix formed of the correlator elements. 
The copy tensor has one input (the physical leg) and 
multiple outputs (two for the case of bond correlators in 
one dimension). This network is equivalent to an MPS, 
where each matrix in the network can be formed from 
the correlators and copy tensors. The copy tensor acts 
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FIG. 12. (a) The copy tensor copies computational basis 
states and subsequently breaks up into disconnected states, 
(b) A generic PEPS in which we expose a single generic rank- 
5 tensor. This tensor network can neither be contracted nor 
sampled exactly and efficiently. However, if the tensor has 
internal structure exploiting the copy tensor then efficient 
sampling becomes possible, (c) The tensor breaks up into a 
vertical and a horizontal rank-3 tensor joined by the COPY- 
dot. Upon sampling computational basis states the result- 
ing contraction reduces to many isolated MPS, each of which 
are exactly contractible, for each row and column of the lat- 
tice. This type of state is known as a string-bond state and 
can be readily generalized (Ref. 14). (d) An even simpler 
case is to break the tensor up into four rank-2 tensors joined 
by a copy tensor forming a nearest-neighbor bond correlator- 
product state, (e) Plaquette correlators, which are outside the 
PEPS class, join up overlapping tensors (in this case rank-4 
ones describing a 2 x 2 plaquette) for each plaquette. Efficient 
sampling is again possible due to the copy tensor. 



by taking the value of the physical leg in the chosen basis 
and copying it to all its legs. As shown in Fig. 11 'b) 
this acts to break up the bonds between the matrices 
while ensuring that neighboring matrices are contracted 
according to the same physical index on a pair of virtual 
legs. The result of this is a product of scalars, shown in 
Fig. 11 'c), which are the correlators for the given input 
configuration used. 

In two dimensions, correlator product states can also 
be considered to be equivalent to a special case of PEPS 
- in this case the copy tensor copies the physical index 
to four different correlator matrices. String bond states 
are another special factorization of PEPS where the con- 
tractible components are themselves MPS, rather than 
a correlator matrix.^ States such as these can then be 
sampled exactly-efflciently, as they consist of product of 
scalar values, each scalar value only depending on a small 
number of physical indices. This is illustrated in Fig. [12 

Since the MPS representation is exhaustive, every cor- 
relator product state has an equivalent matrix product 
state. For certain systems, however, it may be far more 
efficient to use the CPS representation. For example, in 
the MPS representation long-range correlations are me- 
diated by the bonds between neighboring sites. Thus, in 
general, a very large bond-dimension may be needed be- 
tween neighboring tensors to capture the entanglement 
between sites spaced by a large distance. However, using 
CPS it is possible to include a correlator directly be- 
tween distant sites, describing the entanglement with a 
small bond dimension since the bond directly links the 
two sites. Correlator product states can thus describe 



systems that do not obey an area law. 



Appendix B: Numerical methods 

1. Calculating expectation values 

Expectation values of the energy and other operators 
are determined from a CPS by using Monte Carlo impor- 
tance sampling. The energy E is given by 



E = {iP\H\<P) = 



Z ata ,W*(.s')(s>\H\s)W(s) 



(Bl) 



where we assume that the wave function is not normal- 
ized, as is generally the case when using CPS. The energy 
can be expressed as a weighted sum E — J^ s P(s)E(s) 
which is as a sum over configurations of the product of 
the local energy E(s) and the probability P(s), given by 



E(s) 



E 



W*(s') 
W*(s) 



(s'\H\s),P(s) 



\W(s)\ 2 

E.|W(s) 



(B2) 



By replacing the Hamiltonian H with any operator O 
its expectation value (i/j\0\ip) can be similarly computed. 
For example, for the TIM described in section |III[ cal- 
culation of the interaction energy terms cr z a z is 
straightforward since the operators are diagonal: the lo- 
cal energy E zz (s) is given by 



E zz (s) 






(«' 



!.j'] e [i+lJ] _l „[iJ] „[iJ+l] 



), (B3) 



where s^'l = ±1. In a translationally invariant system, 
the sum over all sites can be dispensed with and the 
exchange energy simply taken with respect to the first 
(or any other) spin and its neighbors. 

To calculate the transverse energy terms ga x , we use 
<j x = <j + + <j_ . The only non-zero terms of the energy es- 
timator E(s) = ^2 S , W (J (s'\H\s), are those for which s' 
differs from s by a single spin- flip, so the local transverse 
energy is given by 



E x (s) 



L 

E 

=1,3 = 



W(s) ' 



(B4) 



that the spin on site (i,j) has 
been flipped.^ The absolute magnetization (\cr z \) of the 
ground state can also be calculated. Note that the ex- 
pectation value of the magnetization (a z ) is always be 
zero, since the global Z2 symmetry of the system means 
the ground state forms equal superpositions of configura- 
tions with all spins flipped. However, the absolute mag- 
netization quantifies how well the spins are aligned with 
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one another, and is found by summing the local absolute 
magnetization 



Ms)l = 



L 2 






\i,3\ 



(B5) 



The probability of a given configuration is never ex- 
plicitly calculated using equation (B2). Instead, the con- 
figurations are visited according to importance sampling, 
using the Metropolis Algorithm.^ New configurations s f 
are generated by randomly flipping spins according to the 
calculated acceptance probability |VK(s / )/W r (s)| 2 . There 
are a total of F x N spin flips per bin, where we call F the 
number of sweeps per bin, and N is the number of spins 
(since we visit each spin sequentially in a given sweep). 
As is typical, we also first perform a few warm-up sweeps 
that do not contribute to the estimates to ensure that 
the random walk through configuration space starts in 
an equilibrium position, and also sample the local expec- 
tation value only every s th sweep, where s is the sample 
rate ~ 10, to allow a larger portion of the configuration 
space to be visited and to ensure that the configurations 
that are sampled are independent from one another. 

At no point during the calculation is knowledge of the 
normalization of the wave function needed. The most 
time-consuming step in the calculation is determining the 
correlator fraction W(s')/W(s). Since the configuration 
weight is given by a simple product of numbers, to calcu- 
late this fraction only the correlators that represent the 
sites where the configurations have changed are required. 
If the Hamiltonian is local then there are only a few con- 
figurations s' that have a non-zero matrix element in the 
local energy, so only a few correlator fraction terms need 
to be calculated for each contribution to the energy. 



2. Energy minimization 

The correlator elements that minimize the energy are 
found using a stochastic minimization method which re- 
quires only the sign of the first derivative of the energy 
with respect to the correlator elements.^ Specifically, the 
first derivative is given by 

dE 



oa 



w = 2((Al%E)-(Al%)(E)) J (B6) 



{i} 



.{*> 



where A} { { } (s) is given by 

A S>( s ) 



1 dW(s) 



(B7) 



This is trivial to compute, since W(&) is simply a product 
of the different correlators C. If each correlator is only 
used once A^; } (s) = 1/C}^ }: whereas if the same corre- 
lator is used for multiple sites (e.g. in a translationally 
invariant system) then Ag*{ } (s) = b^/C}^: }: where frw 



is the number of times the correlator Cs {i} appears in the 
product for the correlator amplitude for configuration s. 
Following the procedure in Ref. |33l the expectation 
value of the derivative is calculated in the same way that 
other operators are calculated for F sweeps in a bin, and 
the derivatives for all correlators are calculated simulta- 
neously. After this every correlator is updated according 
to 



C {i} 



C {i} 



rS(k)sign 



dE 



da 






(B8) 



where r is a random number between and 1, and S(k) 
is the step-size for a given iteration k. Unlike the New- 
ton method, the second derivative is not required, con- 
siderably simplifying the calculation. This method has 
been found to give better convergence than the Newton 
method, since it avoids the sizable statistical errors in 
the second derivative and also because it does not ex- 
actly follow the direction of steepest decent which can be 
more efficient. 15 

To achieve convergence, F, the number of sweeps in a 
bin is increased every iteration, and the step size S(k) is 
reduced. For every iteration /c, F is given by F$k, and 
the number of bins is increased as C = Cok per iteration 
to slow the cooling rate, where typically F and Go ~ 10. 
The step size is reduced per iteration as S = Sok -11 . We 
find best results with r] between 0.75 and 0.9. After an 
initial run with a relatively large step size to get close 
to the minimum, the resulting correlators are then used 
as a starting point for a new run of iterations (i.e. k is 
reset) where Fq and Go remain unchanged, but So is de- 
creased. Depending on the system parameters this can 
reduce the energy further, while for other regimes (i.e. 
away from criticality) the energy has already converged 
after the first run. After the minimization is complete, 
the procedure is repeated for a single iteration with zero 
step size and large F and G, to obtain an accurate esti- 
mate of the expectation values. 

For bond correlators with r max < 3 and for 2 x 2 pla- 
quette correlators, we start with a uniform state, and 
perform initial minimization with Fo = 10, Go = 10, 
S = 0.02, r] = 0.9 for 25 iterations. For 3 x 3 and 4 x 4 
plaquette correlators, the initial correlator is built us- 
ing the 2x2 plaquette correlator that results from the 
first round of minimization. For bond correlators with 
r max > 3, we build approximate new correlators using 
the shorter range correlator elements distributed among 
all bond lengths. For all correlator types, the iteration 
counter is then reset, and the minimization procedure 
is performed again using the minimized correlator as a 
starting point, and So = 0.005 for 30-40 iterations (until 
convergence is obtained). Accurate estimations of the ex- 
pectation values are then obtained by performing a final 
calculation with one iteration having a much larger num- 
ber of sweeps per bin (for example F = 10, 000, G = 100, 
6 = 0). A typical minimization routine for plaquette 



correlators is shown in Fig. 13 
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For bond correlators, the best results are found when 
the step size 6 is a function of bond length as well as 
iteration number k. For each bond in the correlator, 
we define the radius r as the maximum displacement 
along one of the lattice indices, and scale the step size 
as <5(r, k) = £(/c)r -0 ' 75 . For the first round of minimiza- 
tion we minimize the bonds for a given radius separately 
in sequence, while in the second round of iterations all 
correlator elements are updated simultaneously. 



FIG. 13. (Color online) Energy during minimization of the 
plaquette correlators for the two-dimensional TIM in a 31 by 
31 system with g = 3.05 (see section III). The initial mini- 
mization routine is performed with 5 — 0.02 for 25 iterations, 
followed by a further round of minimization with 5q = 0.005 
for 30 iterations (although not shown, for 4 x 4 plaquette cor- 
relators the further round of minimization has 40 iterations) 
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